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We describe a new polynomial time quantum algorithm that uses the quan- 
tum fast fourier transform to find eigenvalues and eigenvectors of a Hamiltonian 
operator, and that can be applied in cases (commonly found in ab initio physics 
and chemistry problems) for which all known classical algorithms require expo- 
nential time. Applications of the algorithm to specific problems are considered, 
and we find that classically intractable and interesting problems from atomic 
physics may be solved with between 50 and 100 quantum bits. 
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Long before Shor's ground-breaking algorithm |jj - and the resulting surge of 
interest in quantum computing - Feynman suggested that a quantum computer 
might be useful for simulating other quantum systems 0. This suggestion was 
based upon the observation that quantum systems are described in a Hilbert 
space whose size grows exponentially with the number of particles. Thus a 
collection of only 100 spin 1/2 particles, each of which could be specified by 
only two complex amplitudes were it isolated, requires a total of 2 100 complex 
amplitudes for its state to be specified completely. This exponential explosion 
severely limits our ability to perform true "ab initio" (first principles) calcula- 
tions; since it is obviously not possible to even describe the state of anything 
but the smallest quantum systems, one must resort to various approximation 
techniques to calculate properties of interest. 

Recent work in quantum computation has revealed various techniques for 
simulating physics on a quantum computer |^||ll||^||4||^||5|, and it has been 
demonstrated that this can in fact be accomplished efficiently, as Feynman sup- 
posed. However, while previous work has described a variety of algorithms 
for initializing a quantum computer into a state corresponding to the state of a 
physical system, for time evolving this state on the computer, and for measuring 
properties of the time-evolved state Q |ll[] |2| j| || , there has been comparatively 
little work done on algorithms which calculate static properties of a physical 
system]^. In particular, of all the questions which one might ask about a quan- 
tum system, there is one most frequently asked and for which one would most 
greatly desire an efficient algorithm: What are the energy eigenvalues and eigen- 
states? In this letter, we provide a quantum algorithm that can find eigenvalues 
and eigenvectors of a Hamiltonian operator in cases that occur frequently in 
problems of physical interest. Moreover, the algorithm requires an amount of 
time which scales as a polynomial function of the number of particles and the 
desired accuracy, whereas all known classical algorithms require an exponential 
amount of time. 

The problem to be solved can be precisely stated as follows. Consider the 
time-evolution operator U — e~^ Ht which corresponds to the Hamiltonian H, 
and an approximate eigenvector V a of U (and thus of H) that can be gener- 
ated in quantum polynomial time, i.e., the machine can be placed into a state 
corresponding to V a using a polynomial number of quantum logic operations. 
Call the true eigenvector V and the true eigenvalue A„. If the state V a satisfies 
the property that |(V a |y)| 2 is not exponentially small - that is, the approximate 
eigenvector contains a component of the actual eigenvector that is bounded by 
a polynomial function of the problem size - then V and X v can be found in time 
proportional to 1/ |(V a |l/)| 2 and 1/e, where e is the desired accuracy. 

Intuitively, what the algorithm does is to resolve the guess into its non- 
negligible components and determine the corresponding eigenvalues. If the op- 
erator U (and thus its eigenvectors) is of exponentially large dimension - which 
it typically is - there are no known classical algorithms that can find even the 
eigenvalues in polynomial time. Although the requirement that there exist an 
initial statevector V a with the specified properties may appear to be overly re- 
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strictive, it is frequently (if not usually) possible to obtain such a guess for 
"real" problems using existing classical techniques. For example, in any physi- 
cal system with discrete energy levels that are not exponentially close together 
near the ground state (such as an atom), if it is possible to obtain classically 
any state vector with expected energy merely less than the first excited state 
(by a non-exponentially small amount), then this state vector must contain a 
non-negligible component of the ground state and - although it may not re- 
motely resemble the ground state - could be used as the approximate state V 
to determine the true ground state and ground state energy in polynomial time. 
Finally, if for some problems it is not possible to obtain classically a guess with 
the desired properties, it may often be the case that the state vector V a may be 
generated using a quantum algorithm, such as quantum simulated annealing. 

We will now describe an algorithm which applies to any U that can be 
implemented in quantum polynomial time, whether or not it represents the time 
evolution operator corresponding to a given Hamiltonian. (It was shown in Q 
that the time evolution operator corresponding to any local Hamiltonian can 
be implemented in polynomial time on a quantum computer.) This first part 
of the algorithm was described independently by Cleve et. al. in |Q to find 
eigenvalues (but not eigenvectors) of unitary operators, in that case because the 
eigenvalues of a particular operator can be used to solve the abelian stabilizer 
problem. To begin, we consider a quantum computer consisting of m+l+w 
qubits, where a total of m qubits (to be called the index bits) are used for an 
FFT, a total of 1 qubits describe the Hilbert space in which the operator U acts, 
and w extra working qubits are required for temporary storage. Let M = 2 m . 
The accuracy of the result will grow as 1/M. Assume that the m index qubits 
are initially in the state |0>and that the 1 qubits are initially in the state V 
(hence, the need for V a to be generated in quantum polynomial time). That is, 
the initial state is 



where the w work qubits are assumed to be |0>unless specified otherwise. 
We perform a tt/2 rotation on each of the m index qubits to obtain the state 



Next, one performs a series of quantum logic operations that transform the 
computer into the state 



|*>=|0> \V a > 



(1) 




M-l 





M-l 



(3) 
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This transformation is accomplished by applying the operation U to the second 
set of 1 qubits (which are initially in the state \V a }) j times. It can be imple- 
mented easily by performing a loop (indexed by i) from 1 to M. Using standard 
quantum logic operations, set a flag qubit to the value |l>if and only if i <j 
and perform the operation U conditioned on the value of this flag. Thus only 
those components of the above superposition for which i<j are effected. Finally, 
undo the flag qubit and continue with the next iteration. After M iterations, 

the state above is obtained. 

At this point, it is helpful to rewrite the state in a slightly different manner. 
Label the eigenvectors of U by the states \<pk) and the corresponding eigenvalues 
with Afc. We can then write 



\V a )=J2^\^ k ) (4) 

k 

in which case the state (||) above can be rewritten as 

1 M-l 

i*) = v^E bw'E c *i^> ( 5 ) 



3=0 k 
M-l 



^E c *E \i)^) 3 \^) ( 6 ) 

k j=Q 

If we write Afc as e luJk and exchange the order of the qubits so that the labels 
\(/>k) appear first, the result is seen then most clearly: 



M-l 

i*} = ^=E c ^>E e,; ^'> ( ? ) 

k 3=0 

It is now self-evident that a quantum FFT performed on the m index qubits 
will reveal the phases LUk and thereby the eigenvalues A&. The quantum FFT 
requires only poly(m) operations, whereas the accuracy of the result will scale 
linearly with M or 2 m . Each frequency is seen to occur with amplitude Ck = 
{V a \4>k}] by performing a measurement on the m index qubits, one thus obtains 
each eigenvalue with probability |c/c| 2 . Only a polynomial number of trials is 
therefore required to obtain any eigenvalue for which Ck is not exponentially 
small. If the initial guess \V a ) is close to the desired state (i.e., |< V a \V >| 2 is 
close to 1), then only a few trials may be necessary. 

Moreover, one obtains the eigenvectors as well: once a measurement is made 
and an eigenvalue Xk is determined, the remaining 1 qubits "collapse" into the 
state of the corresponding eigenvector. Of course, the state \4>k) is in some sense 
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"trapped" inside the computer. But since it is impossible to store as classical 
information the 2 l phases associated with the state, one cannot possibly hope 
to do better. However, one is likely to be interested in various properties of 
the eigenvectors, and these can be determined by making various measurements 
on the state. For ab initio quantum calculations, easily obtainable properties 
include those of greatest interest: charge density distributions, correlation func- 
tions, momentum distributions, etc. See for a discussion of how relevant 
physical information can be extracted efficiently from the quantum computer. 

We now consider more precisely how to find the eigenvectors and eigenvalues 
of a "real" Hamiltonian. Generally, one wishes to find energy eigenstates for a 
Hamiltonian of the form 



H = Y J {T i + V i )+Y,Vij (8) 

z—l i>j 

where n is the number of particles, Ti is the kinetic energy, Vi is the external 
potential, and Vij is the interaction between the particles. However, there is no 
reason why these techniques cannot be applied to a different Hamiltonian or to 
one containing additional terms, as long as the Hamiltonian can be separated 
into a sum of local interactions (that is, a sum of terms which act upon only 
k qubits, where k is independent of the number of particles n). (In atomic 
problems, for example, one might include effective interactions such as spin-orbit 
coupling or nuclear finite size effects). Because the Hamiltonian is Hermitian, 
we apply the steps described above to the time evolution operator U (t) = e~ lHt , 
which is unitary and has the same eigenvalues and eigenvectors. This time 
evolution operator is generated using the technique described in ||; the key 
idea is to write 



H = Y J H l (9) 

U(t) = e- lHt = (e- iH ^e- iH2 ^...e- iH ^) m +y"[H i ,H j ]-l— + ... (10) 

i>j 

where each Hi acts on only k qubits at a time. (In the Hamiltonian above, 
each Hi represents one of the terms Tj, Vi, or Vij). Let Ui — e~ tHi ™. Each term 
Ui can be implemented efficiently, because it acts in a space of only k quantum 
bits, where k is small. For large enough m, the second term on the right (and 
the higher order terms) approaches zero. It is therefore possible to generate 
U(t) by acting on the state with each Ui in series, a total of m times. A careful 
analysis || reveals that in order to simulate U (i) with an accuracy e, one needs 
to apply 0(t 2 /e) quantum logic operations.^] 

2 Since U(t) has the same eigenvalues and vectors for all t, this might lead one to falsely 
conclude that the number of operations necessary to find the eigenstates to a given accuracy 
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For a specific problem, the form of the matrices Ui depends greatly on the 
basis set chosen to describe the Hilbert space. Moreover, the choice may strongly 
impact the size of the basis required to describe the system accurately. In the 
usual first quantized representation, each particle is described by a series of 1 
qubits representing a single particle wave function. The system as a whole is 
thus represented with n*l qubits. (It is also possible to use a second quantized 
representation, which may be more efficient for certain problems; see p"l|]). For 
the Hamiltonian above, the matrices Ui can be implemented in a particularly 
efficient manner by using either position space or momentum space for the single 
particle basis, and switching between the two via quantum FFTs. However, for 
most problems, these are not the most effective choices to represent the energy 
eigenstates. Other sets of basis states are generally more efficient and are fre- 
quently employed in conventional classical computations: one possible example 
might be wavelets; another common choice might be single electron solutions for 
an effective potential. As long as the single particle basis is of a fixed size (and 
the reason why we choose a more complicated basis set is for the explicit pur- 
pose of keeping it small) , then the operators Ui can always be calculated in the 
chosen basis and implemented using 0(d 4 ) operations, where d is the dimension 
of the single particle basis set Q. Thus one finds that it is possible to apply 
these quantum algorithms using the more elaborate choices of basis sets that are 
commonly employed in conventional ab initio calculations. (Because there exists 
a fast quantum wavelet transform it may be that a wavelet basis turns out 
to be particularly useful). On the other hand, there is a trade-off between mem- 
ory and speed. By using the position or momentum space representation, one 
needs only 0(poly(k)) = O (poly (log d)) operations to perform each Ui] however 
a large number of qubits are required to describe the eigenstates accurately. By 
choosing a more elaborate basis set, one can vastly reduce the required number 
of qubits, but a much larger number of quantum logic operations 0(k 4 ) may be 
necessary to implement each f/j. Thus one finds that, just as with conventional 
computations, the choice of basis sets in the quantum computation will depend 
upon the specific problem at hand and the specific capabilities of the actual 
computing machine. 

Normally, the initial state V a will be the result of a classical calculation, for 
example, a Hartree-Fock calculation or configuration interaction calculation. 
Any ab initio technique which results in a known wave function can be used. 
(Note that this does not include those techniques which utilize density functional 
theory, as we require a wavefunction, not simply a charge density distribution) . 
If the input wave function is not already symmetrized or antisymmetrizcd, we 
can use the algorithms described in |ll| to do so efficiently. 

Finally, we consider state-of-the-art ab initio calculations of atomic energy 

could be reduced by choosing a shorter length of time t for the operator U(t). However, the 
algorithm requires one to calculate U M , and since U(t) M = U(Mt), one sees that U = U(t) 
must be calculated with greater precision if U M is to be calculated for a fixed precision. In 
fact, since the eigenvectors are determined with a precision proportional to M, the number of 
quantum logic operations required to calculate the energy eigenstates to a precision e is seen 
to scale as e~ 2 . 
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levels in order to compare the quantum algorithm described above with known 
classical techniques. Problems from atomic physics serve as a particularly good 
benchmark because extremely accurate experimental data is widely available. 
The quantum algorithm corresponds most closely to what is known as "complete 
active configuration interaction" or "full configuration interaction" techniques, 
because the many-particle basis set includes all possible products of single par- 
ticle basis vectors. This approach is most valuable in situations where the cor- 
relation energy is large and where many "configurations" are of similar energy 
(this typically occurs when many electrons are in open shells). Unfortunately, 
it is difficult to state precisely the minimum size problem for which the quan- 
tum calculation surpasses the best classical calculations, because a variety of 
sophisticated techniques are used to avoid the exponential explosion in basis 
states. That is, the most accurate classical calculations do not employ directly 
the "full configuration interaction" method. Based on ||l0| , however, we esti- 
mate that a calculation of the energy levels of B (5 electrons), using roughly 
20 angular wavefunctions and 40 radial wavefunctions per particle - for a total 
of 800 single particle wave functions and therefore 800 5 «10 15 full many-body 
basis states - may provide more accurate results than any classical calculation 
performed to date. At the very least, such a calculation would reveal scientifi- 
cally interesting (and classically unobtainable) results with respect to electron 
correlation energies in B and the relative importance of various orders of excited 
configurations. 

A quantum calculation of the B ground state, using a basis set as described 
above, can be accomplished with 60 qubits: f per particle to represent the state 
of the atom (for a total of 50 qubits), 6 or 7 qubits for the FFT, and a few addi- 
tional "scratch" qubits^]. Unfortunately, the two particle operators (generated 
by the coulomb attraction between pairs of electrons) take place in a subspace of 
dimension (2 10 ) 2 ; they therefore are represented by matrices with 2 40 elements. 
Implementing such an operator by brute force is likely to remain intractable for 
the foreseeable future. However, it may be possible to perform the necessary 
transformation using a quantum algorithm. One possible technique is to change 
basis sets: by representing the interacting particles in position space, instead of 
with the orbital basis set, it is easy to calculate the coulomb terms. Thus one can 
transform each particle into position space separately (requiring a small number 
of quantum logic operations), perform the time evolution corresponding to the 
coulomb interaction, and then transform back. Unfortunately, a position space 
representation will require many more qubits. We estimate that 30 qubits per 
particle (10 per dimension, for a real space grid of 1024x1024x1024 per particle) 
will more than suffice. Because these 30 qubits are required only temporarily 

3 The number of qubits required for the FFT is not as large as one might at first suppose, 
based on the earlier statement that the accuracy scales linearly with the size of the FFT. This 
statement is true only for a fixed U. By changing U - in particular, by increasing the length 
of time t in U(t) - one can obtain the eigenvalues to arbitrary precision using a fixed number 
of FFT points. However, the number of points in the FFT must be sufficiently large so as to 
seperate the frequencies corresponding to distinct eigenvectors. This is how the estimate of 6 
or 7 qubits (64 or 128 FFT points) is made. 
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for the 2 particles whose interaction we are considering at any particular stage 
in the algorithm, the new efficient algorithm requires a total of 2 x 30 qubits 
(for the interacting particles), an additional 3 x 10 qubits (for the remaining 
particles), and the same 10 qubits for the FFT and work space. It thus ap- 
pears that in order to realistically perform an "interesting" calculation using 
the algorithms described previously, one will need a quantum computer with 
approximately 100 qubits. Of course, the possibility remains that an efficient 
algorithm for implementing the coulomb interaction could be invented that does 
not require additional working space. 

In conclusion: we have provided a new quantum algorithm which can be used 
to find eigenvectors and eigenvalues of a Hamiltonian operator. The algorithm 
provides an exponential speed increase when compared to the best known classi- 
cal techniques. Problems from atomic physics may be the best place to perform 
the first real calculations, both because accurate experimental data is available 
to verify the resulting calculations, and because the parameters involved appear 
to be within the foreseeable range of small quantum computers. We estimate 
that 50 - 100 qubits arc sufficient to perform "interesting" calculations that are 
classically intractable. Finally, we suggest a couple of interesting questions 
which remain open. First; although we have made estimates regarding numbers 
of required qubits, it would be interesting to calculate accurately the number of 
quantum logic gates required to do an "interesting" problem. Second, a detailed 
analysis of the effects of errors would be worthwhile, as would an analysis of 
error correcting codes in this context. 
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